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The  time  series  data  collected  during  the  HIFiRE-1  flight  experiment  was  sampled 
unevenly  due  to  telemetry  dropouts  and  the  sampling  scheme  selected.  Sampling  unevenly 
results  in  frequency  translation  of  power  to  artificial  sidelobes.  These  sidelobes  distort  the 
real  and  imaginary  components  of  the  Fourier  transform  such  that  the  spectrum  of  the 
sampled  data  no  longer  represents  the  spectrum  of  the  physical  process  generating  the  data. 
These  sidelobes  can  be  eliminated  by  resampling  the  data  onto  an  evenly  spaced  grid  at  the 
cost  of  under  predicting  the  power  of  the  higher  frequency  components  in  the  data.  In  this 
paper,  a  compensation  procedure  is  developed  to  recover  the  power  loss  caused  by 
resampling.  This  compensation  procedure  is  suitable  for  stochastic  time  series  data  having 
red-noise  spectra  such  as  the  pressure  fluctuations  recorded  underneath  laminar  and 
turbulent  boundary  layers. 


Nomenclature 

a j  =  regression  coefficient 

a50  =  overlap  correlation  constant 

b,  =  regression  coefficient 

Cxy  =  coherency 

D  =  location  data  point  of  a  partition 

/  =  frequency,  Hz 
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f0  =  fundamental  frequency,  Hz 

fc  =  Nyquist  frequency,  Hz 

Gxx  =  autospectrum.  Pa2  /  Hz 

Gxy  =  cross-spectrum.  Pa2  /  Hz 

K  =  number  of  partitions 

Keff  =  effective  number  of  partitions 

L  =  number  of  data  samples  in  a  partition 

N  =  number  data  samples  in  a  time  series 

R  =  sum  of  squares 

t  =  time,  sec 

tf  =  time  used  to  account  for  different  time  origins  of  time  series  used  in  cross  spectral  analysis,  sec 

tL  =  time  at  the  end  of  a  partition,  sec 

Tp  =  period,  sec 

Tr  =  sampling  interval,  sec 

vv„  =  Hanning  window  function  coefficients 

y„  =  time  series  data 

er2  =  variance 

r  =  time  shift,  sec 

co  =  angular  frequency,  rad 


I.  Introduction 

THE  Hypersonic  International  Flight  Research  Experimentation  (HIFiRE)  program  is  a  hypersonic  flight  test 
program  executed  by  the  Air  Force  Research  Faboratory  (AFRF)  and  the  Australian  Defence  Science  and 
Technology  Organisation  (DSTO).1' 2  Its  purpose  is  to  develop  and  validate  technologies  critical  to  next  generation 
hypersonic  aerospace  systems.  Candidate  technology  areas  include,  but  are  not  limited  to,  propulsion,  propulsion  - 
airframe  integration,  aerodynamics  and  aerothermodynamics,  high  temperature  materials  and  structures,  thermal 
management  strategies,  guidance,  navigation,  and  control,  sensors,  and  system  components  such  as  munitions, 
submunitions  and  avionics.  The  HIFiRE  program  consists  of  extensive  ground  tests  and  computation  focused  on 
specific  hypersonic  flight  technologies.  Each  technology  program  is  designed  to  culminate  in  a  flight  test. 

The  first  science  flight  of  the  program,  designated  as  HIFiRE-1,  launched  March  22,  2010  at  Woomera 
Prohibited  Area,  Australia.3  This  flight  was  dedicated  to  two  aerothermal  experiments  including  laminar -turbulent 
boundary-layer  transition  on  a  7-degree  half-angle  cone  with  nose  bluntness  of  2.5  mm  and  turbulent  shock- 
boundary-layer  interaction  on  a  33-degree  flare/cylinder  configuration.  Prior  to  flight,  a  series  of  experimental  and 
numerical  studies  were  performed  to  aid  the  design  of  the  aerothermal  experiments  such  as  sensor  selection  and 
placement,  and  vehicle  geometry  and  materials.415 
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The  cone,  cylinder,  and  flare  surfaces  of  the  HIFiRE-1  payload  were  instrumented  with  many  thermocouples, 
heat  transfer  gauges,  and  pressure  transducers.  During  flight,  the  data  from  these  transducers  were  digitized  using  a 
sampling  scheme  with  a  non-constant  time  interval  between  samples.  This  sampling  scheme  was  used  in  order  to 
achieve  maximum  stable  throughput  for  a  variety  of  sensors  on  a  single  flight  computer.  Although  this  sampling 
scheme  maximized  data  acquisition,  special  care  must  be  taken  when  performing  spectral  analysis  on  such 
unevenly-sampled  data.  This  was  especially  true  for  the  HIFiRE-1  high -band  width  pressure  data  that  were  sampled 
at  rates  up  to  60  kHz.  In  addition  to  the  sampling  scheme,  random  data  dropouts  sometimes  occurred  in  the 
telemetry.  Generally,  these  events  were  infrequent,  and  the  dominant  feature  of  data  sampling  was  the  uneven 
sampling  rate.  The  current  paper  summarizes  the  methods  used  to  analyze  the  high  bandwidth  pressure  data,  and 
quantifies  spectral  error  introduced  by  uneven  sampling. 

A  system  consisting  of  a  transducer,  antialiasing  filter,  and  sampling  scheme  is  inherently  nonlinear  even 
though  the  transducers  and  antialiasing  filters  were  linear  devices.  This  nonlinearity  is  a  direct  consequence  of  the 
mutual  dependence  between  frequencies  in  the  frequency  space  of  unevenly  sampled  time  series  data,  and  is 
dependent  on  only  the  sampled  times.  The  spectral  analysis  of  data  measured  by  such  a  system  is  inherently 
complicated  and  requires  special  methods  such  as  least  squares,  "  interpolation,-  slotted  resampling,  ’  -  or 
continuous  time  models.26  -7  Approaches  based  on  these  methods  have  been  extensively  studied  and  developed  in 

17  20  21  28  29  30  31 

other  disciplines  such  as  astronomy,  paleoclimatology,-  seismology,-  biomedical  engineering,-  ’  genetics, 
and  laser  Doppler  velocimetry.24 

The  least  squares  method,  also  referred  to  as  the  Lomb-Scargle  method,  determines  the  Fourier  transform  by 
estimating  the  different  sinusoidal  components  of  a  time  series.  While  more  complicated  to  develop,  the  least 
squares  method  is  robust  requiring  very  little  user  input,  and  can  be  applied  to  any  arbitrary  sampling  scheme 
including  random  sampling.  Moreover,  the  least  squares  approach  has  been  shown  to  be  more  fundamental  than  the 
traditional  Fourier  methods  since,  for  evenly  sampled  time  series  data,  the  two  methods  are  exactly  identical.  For 
unevenly  sampled  time  series  data,  the  two  methods  are  identical  within  the  numerical  error  generated  when 
estimating  the  Fourier  transform  by  numerical  integration.  As  mentioned,  the  spectral  distortion  caused  by  sampling 
unevenly  is  independent  of  the  data  and  depends  only  on  the  sampling  scheme.  Therefore,  estimating  the  power 
spectrum  from  unevenly  sampled  data  results  a  distorted  spectrum  that  does  not  necessarily  represent  the  spectrum 
of  the  physical  process.  The  least  squares  method  is  not  suitable  for  stochastic  time  series  data  with  red-noise 
spectra  such  as  the  surface  pressure  fluctuations  underneath  a  laminar  or  turbulent  boundary  layer  because  the 
analysis  yields  a  distorted  spectrum.  The  method  can  be  useful  for  accurately  computing  the  Fourier  transform  of 
unevenly  sampled,  sinusoidal  periodic,  complex  periodic,  and  almost  periodic  data.  There  are  some  periodic  data 
recorded  during  the  HIFiRE-1  flight  experiment  for  which  the  least  squares  method  is  the  preferred  method  for 
computing  the  Fourier  transform. 

In  general,  the  results  obtained  using  interpolation  methods,  slotted  resampling  methods,  and  continuous  time 
models  introduce  unwanted  distortion,  are  not  robust,  and  can  require  significant  user  input  that  is  dependent  on  the 
nature  of  the  time  series.  Specifically,  interpolating  methods  such  as  nearest  neighbor  or  spline  interpolation  are 
known  to  cause  an  underestimate  of  high  frequency  components  independent  of  the  employed  interpolation 
scheme.-  ’  ’  Moreover,  results  obtained  using  kernel  based  interpolation  methods  are  dependent  on  the 

interpolation  kernel  and  kernel  parameters  selected  such  as  sine,  Gaussian,  Laplacian,  or  rectangular  kernels,  and 
bandwidth,  mainlobe  width,  and  window  width  kernel  parameters.16  Unlike  the  least  squares  method,  these  methods 
are  useful  for  estimating  the  spectral  content  of  unevenly  sampled,  stochastic  signals  having  red-noise  spectra. 

The  use  of  slotted  resampling,  which  estimates  the  autocorrelation  function  by  summing  the  product  pairs  XjXj 
into  bins  lagged  by  f,  -  tp  cannot  be  employed  for  an  arbitrary  sampling  scheme  since  sufficiently  large  gaps  in  the 
data  can  cause  the  method  to  fail.  These  gaps  could  be  filled  in  using  interpolation.  The  interpolation  process,  if 
needed,  and  the  binning  process  both  independently  contribute  unwanted  distortion  to  the  time  series  data. 
Regardless  of  these  negative  attributes,  the  slotted  resampling  method  may  be  a  useful  method  for  analyzing  the 
HIFiRE-1  data.  The  use  of  this  method  for  analyzing  the  HIFiRE-1  data  is  still  under  investigation  and  will  not  be 
discussed  further  in  this  paper. 

The  objective  of  this  paper  is  to  identify  and  outline  a  method  capable  of  accurately  estimating  the  autospectral 
density  function  from  the  unevenly  sampled,  surface  pressure  data  underneath  laminar,  turbulent,  and  shock- 
influenced  boundary  layers  collected  during  the  HIFiRE-1  flight  experiment.  This  objective  requires  some 
background  information,  which  is  developed  prior  to  outlining  an  estimation  procedure  for  the  autospectral  density 
function.  In  Section  II,  sampling  theory  is  developed  for  unevenly  sampled  time  series  data.  Several  useful  tools 
are  identified  from  this  development  including  the  spectral  window  function  and  an  alternate  definition  of  the 
Nyquist  frequency  that  holds  for  unevenly  sampled  time  series  data.  Additionally,  the  nonlinear  nature  of  the 
measurement  and  the  source  of  the  distortion  in  the  Fourier  transform  of  unevenly  sampled  measurements  are 
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discussed.  In  Section  III,  the  discussion  applies  the  methodology  developed  in  Section  II  to  identify  the  nature  of 
the  HIFiRE-1  sampling  scheme.  The  Nyquist  frequencies  for  the  unevenly  sampled,  surface  pressure  data  recorded 
during  the  flight  experiment  are  given  and  compared  to  the  cutoff  frequencies  of  the  antialiasing  filters.  The  range 
in  frequencies  without  distortion,  and  thus  suitable  for  analysis,  is  identified  for  each  high  bandwidth  pressure 
transducer.  In  Section  IV,  the  effects  of  resampling  the  unevenly  spaced  time  series  data  onto  an  evenly  spaced  grid 
via  linear  interpolation  are  discussed.  In  this  section,  it  is  shown  that  resampling  simplifies  the  process  for 
estimating  the  autospectral  density  function  of  unevenly  sampled  data  by  replacing  the  distortion  caused  by  the 
nonlinear  sidelobes  with  an  underestimate  in  power  of  the  higher  frequencies  in  the  autospectral  density  function.  In 
Section  V,  a  compensation  procedure  is  outlined  for  recovering  the  power  loss  caused  by  resampling.  This 
compensation  method  is  suitable  for  estimating  the  autospectral  density  function  of  stochastic  time  series  data 
having  red-noise  spectra  such  as  the  surface  pressure  fluctuations  recorded  underneath  laminar  and  turbulent 
boundary  layers  during  the  HIFiRE-1  flight  experiment. 


II.  Sampling  Theory 


A  block  diagram  representing  the  sampling  process  of  a  continuous  function  is  given  in  Fig.  1.  In  Fig.  1,  the 
signal  x(t)  represents  the  output  from  a  measurement  device  such  as  a  pressure  transducer  and  is  a  continuous 
function,  x'(t )  is  a  band  limited  function  assumed  to  have  the  same  amplitude  and  phase  spectrum  as  x(t)  for 
frequencies  below  the  cutoff  frequency  of  the  antialiasing  filter,  s(t)  is  the  sampling  function,  and  y(tn)  is  a  set  of 
discrete  samples  of  x'(t)  sampled  at  times  tn.  The  function  x(t)  is  filtered  before  sampling  so  that  the  frequency 
content  of  x'(t)  is  limited  to  frequencies  below  the  Nyquist  frequency,  defined  as 


fc  = 


1 

2A  t 


(1) 


where  At  represents  the  sampling  interval.  According  to  the  sampling  theorem,  x'(t )  can  be  reconstructed  from 
y(tn)  without  information  loss  or  distortion  if  the  sampling  rate  is  greater  than  or  equal  to  twice  the  Nyquist 
frequency.  If  the  sampling  rate  is  less  than  twice  the  Nyquist  frequency,  spectral  leakage  from  frequencies  greater 
than  the  Nyquist  frequency  wrap-around  and  combine  with  the  frequency  content  of  frequencies  below  the  Nyquist 
frequency.  This  problem  is  known  as  aliasing  and  in  general  cannot  be  removed. 

The  sampling  function  is  defined  as 


s(t)  =  cln=$nSlt  tn) 

Ln=lwn 


(2) 


where  S(t  —  tn)  is  the  Dirac  delta  function,  t„  are  the  sampled  times,  N  is  the  number  of  data  points,  and  C  and  w„ 
are  constants.  The  discrete  signal  y(tn)  is  given  as 

y(tn)  =  s(t)x\t)  =  ~2n=i  x' (tn)8(t  -  tn)  (3) 

where  the  constants  C  and  wn  were  set  equal  to  1 .  The  discrete  Fourier  transform  of  Eq.  3  is  given  as 

Yin  =  sin  *  x'io  (4) 

where  the  Fourier  transform  is  defined  as 


Y(fk)  =  C  yit)eWdt  «  yn  e2nifktnht.  (5) 

Typically,  the  sampling  interval  is  moved  outside  the  summation  because  it  is  constant  for  an  evenly  sampled  time 
series.  In  this  paper,  the  sampling  interval  is  not  necessarily  constant  and  the  sampling  interval  is  left  within  the 
summation.  In  Eq.  4,  Y(f)  is  the  Fourier  transform  of  y(tn),  X’ if)  is  the  Fourier  transform  of  x'(tn),  and  S(f)  is 
the  Fourier  transform  of  s(t)  and  is  referred  to  as  the  normalized  spectral  window  function.  According  to  Eq.  4,  and 
as  pointed  out  by  Deeming32,  "the  pathology  of  the  data  distribution  is  all  contained  in  the  spectral  window,  which 
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can  be  calculated  from  the  data  spacing  alone,  and  does  not  depend  directly  on  the  data  themselves.”  The  discussion 
will  now  focus  on  identifying  the  source  of  the  distortion  in  the  Fourier  transform  of  unevenly  sampled  time  series 
data  using  the  normalized  spectral  window  function. 


*■  x'(‘)  — — ♦ao 


s(t) 

Figure  1.  Block  diagram  depicting  the  sampling  process  of  a  continuous  function. 


The  power  spectrum  of  the  spectral  window  function,  defined  as  |S(/)|2,  is  shown  in  Fig.  2  for  an  evenly 
sampled  time  series  with  a  sample  rate  of  100  Hz.  From  Fig.  2,  the  normalized  spectral  window  function  has  peaks 
with  amplitude  of  1  centered  at  2 nfc  where  n  =  0,  1,  2  ...  co.  Hence,  the  normalized  spectral  window  function  is  a 
periodic  function  with  period  2 fc.  This  periodicity  is  independent  of  the  sample  rate. 

The  power  spectrum  of  the  normalized  spectral  window  function  for  an  unevenly  sampled  time  series  sampled  at 
a  sampling  rate  of  200  Hz  but  with  gaps  in  the  data  is  shown  in  Fig.  3a.  The  power  spectrum  in  Fig.  3a  is  similar  to 
the  power  spectrum  in  Fig.  2  with  peaks  centered  at  2 nfc,  however,  there  are  additional  peaks  centered  at  2,  6,  and  10 
Hz,  as  shown  in  Fig.  3b.  These  additional  peaks  are  commonly  referred  to  as  sidelobes.  The  magnitude  and 
location  of  the  sidelobes  are  dependent  on  the  sampling  scheme.  As  will  be  shown,  these  sidelobes  cause  the  real 
and  imaginary  components  of  Y  (/)  and  X'(f)  to  be  different.  This  difference  is  referred  to  within  this  paper  as  the 
distortion  in  Y (/)  caused  by  the  sampling  scheme. 

The  period  of  the  spectral  window  function  for  an  unevenly  sampled  time  series  is  approximately  2 fc,  as  shown 
in  Fig.  3a.  The  spectral  window  function  of  an  unevenly  sampled  time  series  is  not  always  exactly  periodic,  as  the 
peaks  centered  at  2 nfc  for  n  >  0  are  sometimes  split  into  multiple  peaks.  The  separation  in  frequency  of  the  split 
peaks  is  very  small  such  that  the  peaks  are  centered  at  approximately  2 nfc.  The  center  of  the  peak(s)  at  2 nfc  is  used 
to  determine  the  Nyquist  frequency  of  an  unevenly  sampled  time  series  since  Eq.  1  is  not  valid  when  the  sampling 
interval  is  not  constant33.  The  meaning  of  the  Nyquist  frequency  determined  from  the  spectral  window  function  is 
unchanged  i.e.  signals  with  non-zero  amplitudes  at  frequencies  exceeding  the  Nyquist  frequency  wrap-around 
polluting  the  Fourier  transform  at  frequencies  below  the  Nyquist  frequency. 

There  is  an  additional  unwanted  artifact  caused  by  aliasing  that  is  inherent  only  to  unevenly  sampled  time  series 
data.  For  example,  let  x'(t)  =  cos(2ntf )  with/=  96  Hz  and  let  y(tn)  be  the  set  of  discrete  samples  of  x'(t )  where 
the  sampling  scheme  has  the  normalized  spectral  window  function  depicted  in  Fig.  3a.  The  Nyquist  frequency  for 
this  sampling  scheme  is  100  Hz.  The  sidelobes  centered  at  2,  6,  and  10  Hz  in  Fig.  3b  are  centered  at  98,  104,  and 
106  Hz  after  convolving  S(/)  with  X’  (f).  The  sidelobes  centered  at  104  and  106  Hz  are  greater  than  the  Nyquist 
frequency  and  wrap-around  to  96  and  94  Hz,  respectively.  Notice  the  peak  in  the  power  spectrum  corresponding  to 
the  input  signal  has  a  frequency  below  the  Nyquist  frequency  and  does  not  wrap-around  but  combines  with  the 
sidelobe  aliased  to  96  Hz.  The  real  and  imaginary  components  of  Y (/)  at  96  Hz  are  distorted  due  to  the  aliased 
sidelobe. 

Another  property  of  interest  is  now  given.  Assume  that  x(t)  is  the  output  of  a  constant  parameter  linear  system. 
If  the  output  of  such  a  system  is  sampled  unevenly,  then  the  overall  system  is  nonlinear.  This  property  is 
demonstrated  with  the  following  example.  Let  x'(t )  =  cos(2ntf )  with /=  23  Hz  and  let  y(tn)  be  the  set  of  discrete 
samples  of  x’(t )  where  the  sampling  scheme  has  the  normalized  spectral  window  function  depicted  in  Fig.  3a.  The 
Fourier  transform  of  x'(t )  will  be  zero  everywhere  except  at  /=  23  Hz.  The  Fourier  transform  of  Y  (/)  according  to 
Eq.  4  is  S(f)  *  X’(f).  The  power  spectrum  of  Y (/)  is  shown  in  Fig.  4.  In  Fig.  4,  the  Fourier  transform  of  Y (/)  is 
computed  by  convolving  S(f)  with  X'(f  )  and  by  computing  the  Fourier  transform  of  y(tn).  The  Fourier  transform 
of  y(tn)  was  computed  using  the  least  squares  Fourier  transform  to  eliminate  error  caused  by  numerical  integration. 
The  least  squares  Fourier  transform  is  shown  to  be  identical  to  the  traditional  Fourier  transform  defined  in  Eq.  5  in 
Appendix  A.  In  Fig.  4,  the  two  methods  for  computing  Y (/)  are  identical.  More  importantly,  the  sampling  scheme 
results  in  the  leakage  of  power  from  23  Hz  to  sidelobes  centered  at  13,  17,  21,  25,  29,  and  33  Hz.  That  is,  a  system 
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sampled  unevenly  causes  frequency  translation  and  is  therefore  nonlinear  regardless  of  whether  the  transducer  has  a 
linear  response. 

Due  to  the  inherent  nonlinearity  caused  by  sampling  unevenly,  the  reconstruction  of  x'(t )  from  s(t)  and  y(tn) 
does  not  in  general  exist.  However,  it  has  been  shown  that  it  is  possible  to  reconstruct  x’(t )  exactly  from  s(t)  and 
y(tn)  if  the  power  spectrum  of  x’(t)  is  sparse,  and  if  the  frequency  content  of  x’(t)  is  below  the  Nyquist  frequency 
evaluated  at  the  mean  sampling  interval34.  This  reconstruction  process  requires  nonlinear  methods  such  as  an 
iterative  deconvolution  and  is  important  for  the  research  area  of  compressed  sensing.  Unfortunately,  the  power 
spectra  of  surface  pressure  fluctuations  underneath  laminar,  turbulent,  or  shock-influenced  boundary  layers  show  a 
continuous  decrease  of  spectral  amplitude  with  increasing  frequency  (red-noise  spectrum)  and  hence  are  not  sparse. 
Therefore,  methods  that  estimate  x\t)  from  s(t)  and/or  y(tn)  have  to  be  used. 


Frequency,  Hz 

Figure  2.  Power  spectrum  of  the  normalized  spectral  window  function  for  an  evenly  sampled  time  series  with  a 
sample  rate  of  100  Hz. 


Figure  3.  (a)  Power  spectra  for  the  spectral  window  function  of  an  unevenly  sampled  time  series,  (b)  Rescaled  view 
of  Fig.  3(a)  showing  the  sidelobes.  The  sampling  interval  was  constant  at  0.005  seconds  with  gaps  between  0.250  to 
0.500,  0.750  to  1 .000,  and  1 .250  to  1 .500  seconds. 
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Frequency,  Hz 

Figure  4.  Power  spectrum  computed  by  convolving  a  sinusoidal  input  at  a  single  frequency  of  23  Hz  with  the 
normalized  spectral  window  function,  and  by  computing  the  Fourier  transform  of  the  unevenly  sampled  time  series. 

The  discussion  has  so  far  focused  on  periodic  data.  The  discussion  is  now  extended  to  stochastic  time  series 
data  such  as  surface  pressure  data  recorded  underneath  laminar,  transitional,  turbulent,  or  shock -impinged  boundary 
layers.  The  power  spectra  of  this  type  of  data  will  have  nonzero  power  at  each  resolvable  frequency.  The 
frequencies  with  nonzero  power  each  contribute  sidelobes  resulting  in  significant  distortion  of  the  Fourier  transform 
due  to  the  sum  contribution  of  so  many  sidelobes.  The  distortion  in  the  Fourier  transform  of  surface  pressure 
measurements  underneath  a  transitional  boundary  layer  is  shown  in  Fig.  5.  Evenly  sampled  surface  pressure  data 
was  acquired  at  a  sample  rate  of  500  kHz13.  The  corresponding  Nyquist  frequency  was  250  kHz.  Two  unevenly 
sampled  time  series  were  generated  from  the  evenly  sampled  time  series  by  deleting  every  sixth  and  every  third  data 
point.  Each  time  series  was  partitioned  into  400  segments.  The  Fourier  transform  of  each  segment  was  computed 
using  the  least  squares  Fourier  transform.  The  Fourier  transform  of  each  time  series  was  determined  by  averaging 
the  Fourier  transform  from  each  segment.  The  power  spectrums  of  the  evenly  sampled  and  two  unevenly  sampled 
time  series  are  shown  in  Fig.  5a.  In  Fig.  5a,  the  distortion  caused  by  the  sidelobes  is  apparent  by  comparing  the 
amplitude  of  the  evenly  sampled  time  series  (ESTS)  with  the  amplitudes  of  the  two  unevenly  sampled  time  series 
(USTS).  The  phase  was  determined  for  each  time  series  by  computing  the  arctangent  of  the  ratio  of  the  imaginary 
and  real  component  of  the  Fourier  transform.  A  comparison  between  the  phase  of  the  evenly  sampled  time  series 
and  one  of  the  unevenly  sampled  time  series  is  shown  in  Fig.  5b.  The  phase  spectrum  from  the  other  unevenly 
sampled  time  series  was  not  included  for  clarity  purposes.  In  Fig.  5b,  the  phases  of  the  evenly  and  unevenly 
sampled  time  series  are  different  due  to  the  distortion  caused  by  the  sidelobes. 


Figure  5.  (a)  Power  spectrum  of  an  evenly  sampled  time  series  (ESTS),  and  two  unevenly  sampled  time  series 
generated  by  deleting  every  sixth  data  point  (USTS  1)  and  every  third  data  point  (USTS  2)  of  the  ESTS,  (b)  A 
comparison  between  the  phase  spectrum  of  the  ESTS  and  USTS.  The  time  series  data  correspond  to  surface 
pressure  fluctuations  recorded  underneath  a  transitional  boundary  layer. 
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III.  The  HIFiRE-1  Sampling  Scheme 


The  sampling  scheme  for  the  high  bandwidth  pressure  transducers  on  the  cone,  flare,  and  cylinder  are  shown  in 
Figs.  6a,  6b,  and  6c,  respectively.  In  this  paper,  the  normalized  time  interval  is  defined  as  At/min(At)  and 
represents  the  duration  between  consecutive  measurements  in  units  of  the  time  interval.  As  shown  in  Fig.  6,  the 
sampling  scheme  used  for  each  high  bandwidth  pressure  transducer  is  periodic.  The  high  bandwidth  pressure 
transducer  designated  as  transducer  12  was  sampled  at  a  lower  rate  than  the  rest  of  the  transducers  located  on  the 
cylinder  to  satisfy  telemetry  bandwidth  constraints,  as  shown  in  Fig.  6d.  As  will  be  shown,  the  sample  rate  for 
transducer  12  was  well  below  the  cutoff  frequency  of  the  antialiasing  filter.  The  sampling  scheme  used  for  the 
transducers  distributed  on  the  flare  consisted  of  48  consecutively  measured  data  points  sampled  evenly  for  1 .2407  x 
10"3  seconds,  as  shown  in  Fig.  6b.  The  fundamental  frequency  for  these  evenly  sampled  blocks  of  data  is  -806  Hz 
and  hence  the  frequency  resolution  is  poor.  However,  the  Fourier  transform  of  these  blocks  of  data  are  free  of  the 
distortion  caused  by  sidelobes. 


Figure  6.  The  sampling  scheme  for  the  high  bandwidth  pressure  transducers  mounted  on  (a)  the  cone,  (b)  the  flare, 
(c)  the  cylinder  (transducers  7,  8,  9,  10,  11,  and  19),  and  (d)  the  cylinder  (transducer  12). 


The  power  spectrum  of  the  normalized  spectral  window  functions  for  the  high  bandwidth  pressure  transducers 
mounted  on  the  cone,  flare,  and  cylinder  are  shown  in  Figs.  7a,  7b,  and  7c,  respectively.  The  power  spectrum  of  the 
normalized  spectral  window  function  is  riddled  with  many  sidelobes.  As  demonstrated  in  the  previous  section,  these 
sidelobes  cause  significant  distortion  of  the  Fourier  transform.  The  mean  Nyquist  frequency,  computed  using  the 
mean  sampling  interval  and  Eq.  1,  are  given  as  17.81,  15.87,  and  22.67  kHz  for  the  transducers  mounted  on  the 
cone,  flare,  and  cylinder,  respectively.  The  cone  transducers  were  bandpass-filtered  at  100  Hz  -  30  kHz,  and  those 
on  the  cylinder  and  flare  were  filtered  at  100  Hz  -  20  kHz.  Signals  on  the  cone  and  flare  were  thus  subject  to  some 
aliasing  at  higher  frequencies,  but  the  cylinder  transducers  should  be  free  of  aliasing. 
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Figure  7.  Power  spectrum  of  the  spectral  window  function  for  the  high  bandwidth  pressure  transducers  mounted  on 
(a)  the  cone,  (b)  the  flare,  (c)  the  cylinder  (transducers  7,  8,  9,  10,  11,  and  19),  and  (d)  the  cylinder  (transducer  12). 


IV.  Linear  Interpolation 


Resampling  irregularly  spaced  data  onto  a  uniformly  spaced  grid  by  interpolation  eliminates  the  nonlinear 
sidelobes,  as  can  be  shown  by  computing  the  spectral  window  function.  However,  interpolation  acts  as  a  low-pass 
filter,  significantly  reducing  the  power  of  the  higher  resolvable  frequencies  of  the  data.  The  amount  of  power  loss 
depends  on  the  interpolation  scheme.  In  general,  every  interpolation  scheme  acts  as  a  low-pass  filter  causing  some 
power  loss36. 

The  grid  spacing  for  an  evenly  spaced  grid,  which  is  commonly  referred  to  as  the  resampling  time,  is  given  as 


where  the  Nyquist  frequency  is  determined  from  the  period  of  the  spectral  window  function.  The  uniform  grid  is 
given  as  tn  =  t0  +  (n  —  l)tr  for  n  =  1,  2,  3....  (tj-  —  t0)/tr  where  t0  and  tf  are  the  initial  and  final  times  of  the 
time  series. 

The  discussion  will  now  focus  on  estimating  the  power  loss  caused  by  resampling  the  HIFiRE-1  data  onto  a  grid 
with  grid  spacing  tr.  The  objective  will  be  to  determine  a  time-series  with  analytically  described  spectral 
characteristics  similar  to  the  expected  HIFiRE-1  pressure  signals.  The  spectral  distortion  caused  by  uneven 
sampling  can  then  be  quantified,  and  this  known  distortion  can  then  be  used  to  determine  the  transfer  function  of  the 
resampling  process. 

Consider  Langevin’s  equation  describing  the  velocity  x(t„)  of  a  Brownian  particle  given  as 
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dx(t )  =  r  x(t)dt  +  d.f(t) 


(7) 


where  r  <  0  and  f  (t)  represents  a  random  process  consistent  with  Brownian  motion.  According  to  Robinson37,  the 
velocity  samples  generated  by  Eq.  7  are  also  generated  by 

x(tn)  =  exp  (-  x(tn_,)  +  £(tn)  (8) 


where 


£(tn)~JV  {o,  1  -  exp  (-2  ^±)}. 


(9) 


In  Eq.  9,  £  represents  the  random  component  of  the  time  series  and  is  assumed  to  have  a  normal  distribution  with  a 
mean  and  variance  of  zero  and  1  —  exp  (—2  tn  respectively.  The  autospectral  density  function  can  be 

determined  by  computing  the  Fourier  transform  of  the  autocorrelation  function  of  Eq.  8  and  is  given  as 


Gxx(D  ~  Go' 


i -pz 


-2 P  COS 


where  G0  is  the  average  spectral  amplitude  and  p  is  given  as 

P  =  exp  (-7). 


(10) 


(11) 


In  Eq.  11,  r  is  a  time  constant  and  represents  the  decay  period  of  the  autocorrelation  function,  and  At  is  the  average 
sampling  interval  given  as 


At  =  ^1.  (12) 

According  to  Robinson37,  Eq.  8  may  be  thought  of  as  a  first  order  autoregressive  model  with  time-varying 
coefficients  and  heteroscedastic  errors. 

Equations  8  and  10  provide  a  way  to  compute  time  series  data  with  a  known  autospectral  density  function  for 
any  arbitrary  sampling  scheme.  Equations  8  through  12  were  validated  by  generating  time  series  data  for  uniformly 
spaced  samples  and  then  comparing  the  autospectral  density  functions  computed  using  Eq.  10  and  computed  directly 
from  the  time  series  data  using  the  procedures  given  in  Bendat  and  Piersol38.  The  comparison  is  nearly  identical,  as 
shown  in  Fig.  8a. 

The  power  loss  caused  by  resampling  the  HIFiRE-1  data  was  estimated  by  generating  time  series  data  using  Eq. 
8  for  the  HIFiRE-1  sampling  scheme.  This  time  series  was  interpolated  onto  an  evenly  spaced  grid.  The 
autospectral  density  function  was  determined  directly  from  Eq.  10  and  compared  with  the  autospectral  density 
function  computed  from  the  interpolated  time  series.  Comparisons  are  shown  in  Figs.  8b,  8c,  and  8d  for  the 
sampling  schemes  used  by  the  pressure  transducers  located  on  the  cone,  flare,  and  cylinder,  respectively.  The  results 
shown  in  Figs.  8b,  8c,  and  8d  are  for  a  linear  interpolation  scheme.  Other  interpolation  methods  were  investigated 
and  yielded  similar  results.  As  shown  in  the  Fig.  8,  the  power  loss  was  greatest  for  the  transducers  located  on  the 
cone,  whereas,  the  transducers  located  on  the  flare  experienced  the  least  power  loss.  It  is  reasonable  that  the 
transducers  on  the  cone  experienced  a  greater  power  reduction  than  the  transducers  on  the  flare  because  the  data 
collected  on  the  cone  were  missing  more  samples  than  the  data  collected  on  the  flare,  as  shown  in  Fig.  6.  Similarly, 
blocks  of  data  riddled  with  many  missing  data  points  such  as  near  the  beginning  of  the  first  and  second  stage  rocket 
burn  will  experience  the  most  power  reduction.  The  uncertainty  for  these  blocks  of  data  will  be  greater  than  blocks 
of  data  with  fewer  missing  samples. 
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Frequency,  kHz 

Figure  8.  (a)  Equation  8  and  Eq.  10  are  compared  for  a  uniformly  spaced  sampling  scheme,  (b)  The  sampling 
scheme  of  the  pressure  transducers  located  on  the  cone  are  used  to  generate  a  time  series  with  Eq.  8.  The  time  series 
is  interpolated  and  the  autospectral  density  function  is  computed  and  compared  to  the  theoretical  autospectral 
density  function,  (c)  The  theoretical  autospectral  density  function  is  compared  to  the  autospectral  density  function 
computed  from  interpolated  data  for  the  sampling  scheme  used  by  the  transducers  mounted  on  the  flare,  (d)  The 
theoretical  autospectral  density  function  is  compared  to  the  autospectral  density  function  computed  from 
interpolated  data  for  the  sampling  scheme  used  by  the  transducers  mounted  on  the  cylinder. 


V.  Compensation  Method  for  Estimating  the  Autospectral  Density  Function 

As  discussed  in  the  previous  section,  interpolating  unevenly  sampled  data  onto  an  evenly  spaced  grid  eliminates 
the  nonlinear  sidelobes  but  acts  as  a  low-pass  filter  reducing  the  power  of  the  higher  frequency  content  of  the  time 
series  data.  In  this  section,  a  method  is  developed  to  compensate  the  reduction  in  power  caused  by  linear 
interpolation.  The  compensation  method  developed  in  this  section  is  very  similar  to  the  method  given  by  Schulz  and 
Mudelsee39.  After  compensation,  the  autospectral  density  function  is  similar  to  the  autospectral  density  function  of 
x\t). 

It  is  assumed  that  the  time  series  data  does  not  have  sinusoidal  periodic,  complex  periodic,  or  almost-periodic 
components  and  is  a  stochastic  time  series  with  a  red-noise  “like”  autospectral  density  function.  Sinusoidal 
components  should  be  removed  from  the  time  series  before  compensation.  The  steps  for  the  compensation  method 
are  as  follows: 

1.  Determine  r  by  fitting  the  measured  data  with  Eq.  8.  There  are  two  equivalent  approaches  for  determining  r  in 
Eq.  8  from  the  measured  data.  The  first  approach  is  to  use  the  estimation  procedure  given  by  Mudelsee40,  which 
utilizes  Brent’s  method.  The  second  approach  consists  of  more  steps  but  utilizes  the  equations  given  in  this  paper. 
To  begin,  interpolate  the  measured  data  onto  a  uniformly  spaced  grid  and  determine  the  autoregressive  coefficient 
given  as 
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(13) 


<Pi  =  exP  (-  7) 

by  fitting  the  interpolated  data  using  least  squares  given  as 


j?n=n+l  ^m-l 


jm=n+l  |  \V 0 

?M 


2]  [<? 

-  u 


Zm=n+1  ^m^m-1 
VM  r  r 
£-im=n+l  A'mA'm—2 


lSm=n+l  %m— lAm— 2  X?n=n+1  -*771—2  1  L^lJ 

In  Eq.  14,  M  is  the  total  number  of  autoregressive  coefficients.  The  time  constant  is  given  as 

At 


Tfit  ■ 


logoi' 


(14) 


(15) 


The  time  constant  given  in  Eq.  15  does  not  actually  correspond  to  r  of  the  measured  data.  Instead,  an  iterative 
procedure  is  required  to  determine  the  appropriate  value  for  r.  To  do  so,  r  is  assigned  an  initial  value  designated  as 
t'.  Next,  compute  a  time  series  using  Eq.  8  for  the  same  sampling  scheme  as  the  measured  data  for  r  =  r'. 
Interpolate  this  time  series  onto  the  same  uniformly  spaced  grid  as  the  measured  time  series  data.  Determine  the 
time  constant  designated  as  f  for  this  time  series  using  Eqs.  14  and  15.  Adjust  r'  so  that  |f  —  Tfit\  — >  0.  The  time 
constant  of  the  measured  data  is  then  estimated  as  r'. 


2.  Compute  an  autoregressive  time  series  using  Eq.  8  for  r  =  r'. 

3.  Interpolate  the  time  series  data  from  Step  2  and  compute  the  autospectral  density  function.  The  autospectral 
density  function  from  this  step  is  designated  as  G’ AR . 

4.  Compute  the  autospectral  density  function  of  the  time  series  data  from  Step  2  using  Eq.  10.  The  autospectral 
density  function  from  this  step  is  designated  as  Gar- 

5.  Compute  the  compensation  factor  given  as 


6.  Determine  the  compensated  autospectral  density  function  given  as 

Gxx(f)  =  WWXX<J )  (17) 

where  G'xx(f )  is  the  autospectral  density  function  of  the  measured  data  interpolated  onto  an  evenly  spaced  grid. 

The  compensation  method  was  checked  in  two  ways.  In  one  approach,  short  blocks  of  flight  data  with  even 
sampling  rates  provided  spectra  undistorted  by  uneven  sampling  time  intervals.  In  the  second  approach,  a  large - 
eddy  simulation  (LES)  of  a  turbulent  shock  boundary-layer  interaction  provided  a  time-series  that  could  be  evenly 
and  unevenly  sampled.  The  evenly-sampled  LES  time  series  provided  a  truth-model  against  which  to  compare  the 
unevenly-sampled,  compensated  spectrum. 

As  shown  in  Fig.  6b,  the  flight  sampling  scheme  on  the  flare  featured  evenly  sampled  blocks  of  data  separated  by 
gaps  with  duration  of  approximately  10.5  sampling  intervals.  The  autospectrum  was  determined  by  averaging  the 
autospectrum  computed  for  1200  evenly  sampled  blocks  of  data.  This  spectrum  is  treated  as  a  known  spectrum  with 
frequency  resolution  of  approximately  800  Hz.  This  known  spectrum  is  designated  as  G(ESTS).  The  autospectrum 
of  unevenly  sampled  data  was  determined  from  the  flare  data  by  partitioning  the  data  into  blocks  with  the  gaps.  As 
shown  in  Fig.  9,  G'  xx(f)  is  significantly  lower  than  G(ESTS)  for  frequencies  greater  than  2  kHz.  A  comparison 
between  Gxx(/)  and  G(ESTS)  shows  that  the  compensation  method  recovers  the  correct  roll  off  and  power  levels. 
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Figure  9.  Comparison  between  G'xx(/ )  and  Gxx(/)  with  a  baseline  spectrum,  designated  as  GUESTS')  for  surface 
pressure  measurements  recorded  by  transducer  21  located  on  the  flare.  The  baseline  spectrum  was  determined  from 
the  evenly  sampled  blocks  of  data  shown  in  Fig.  6b. 

The  compensation  method  was  also  validated  using  numerical  time  series  data.41  This  numerical  time  series  was 
sampled  evenly  at  a  sample  rate  of  192.9  MHz.  The  numerical  time  series  was  resampled  to  match  the  sampling 
scheme  used  for  the  transducers  located  on  the  cone.  The  numerical  time  series  was  digitally  low-pass  filtered  prior 
to  resampling  in  the  same  pass-band  as  the  HIFiRE-1  data.  Resampling  was  accomplished  by  sampling  the  time 
series  using  the  sampling  scheme  shown  in  Fig.  6a.  Since  the  numerical  time-series  was  sampled  at  discrete  time 
intervals,  the  uneven  sampling  rate  was  similar  to  the  HIFiRE  sampling  rate,  but  did  not  exactly  reproduce  it.  The 
spectral  window  functions  for  the  sampling  scheme  of  the  cone  transducers  and  the  resampled  numerical  time  series 
are  shown  in  Fig.  10.  As  shown,  the  spectral  window  functions  are  very  similar,  thus  the  distortion  caused  by 
sidelobes  in  the  autospectral  density  function  should  be  very  similar  for  the  two  time  series.  In  Fig.  11,  the 
autospectral  density  function  was  determined  for  the  evenly  spaced  numerical  time  series,  the  unevenly  spaced 
numerical  time  series  interpolated  onto  an  evenly  spaced  grid,  and  the  unevenly-spaced  time  series,  interpolated  and 
compensated.  As  shown  in  Fig.  11,  the  compensation  method  recovers  the  roll-off  and  power  loss  caused  by 
interpolating  the  unevenly  sampled  numerical  time  series  onto  an  evenly  spaced  grid.  This  numerical  experiment 
and  the  comparison  of  unevenly  to  evenly-spaced  flight  data  confirm  the  validity  of  the  compensation  method. 
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Figure  10.  Comparison  between  the  spectral  window  functions  of  the  sampling  scheme  used  in  flight  by 
the  transducers  located  on  the  cone  with  the  resampled  numerical  time  series.  The  frequencies  were 
normalized  by  the  Nyquist  frequency  of  each  sampling  scheme. 
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Figure  11.  Comparison  between  the  autospectral  density  function  of  the  compensated  spectrum,  the 
evenly  sampled  numerical  time  series,  and  the  resampled  numerical  time  series  interpolated  onto  an 
evenly  spaced  grid. 

A  comparison  between  the  dimensional  and  nondimensional  G'xx(/)  and  Gxx(/)  for  surface  pressure 
measurements  recorded  at  different  times  during  the  ascent  by  transducers  7  and  9  located  on  the  cylinder  are  shown 
in  Figs.  12  and  13.  These  transducers  were  located  respectively  8  and  6  cm  upstream  of  the  flare.  The  frequency 
and  autospectral  density  function  were  normalized  using  the  99  percent  boundary-layer  thickness,  and  the  dynamic 
pressure.  These  normalizing  quantities  were  provided  by  Robert  Yentsch  of  The  Ohio  State  University.  The  flight 
conditions  for  the  times  shown  in  Figs.  12  and  13  are  summarized  in  the  table  below.  Times  prior  to  t=6  seconds 
correspond  to  the  first-stage  boost  phase.  Times  between  t=6  and  t=15  seconds  correspond  to  the  coast  phase 
between  first-stage  burnout  and  second-stage  ignition. 

The  flight  data  show  a  rolloff  in  amplitude  with  increasing  frequency,  typical  of  that  seen  in  shock-boundary 
interactions.42  As  shown  in  the  figure,  the  compensation  method  does  not  noticeably  alter  power  levels  for 
frequencies  below  1.5-3  kHz  but  becomes  more  prominent  as  the  frequency  increases.  The  compensated  spectra 
roll  off  at  a  rate  of  about  8  dB/octave.  This  rolloff  rate  is  slightly  lower  than  that  observed  in  Ref  41,  but  similar  to 
other  data  cited  in  the  same  reference.  These  data,  combined  with  checks  against  numerical  experiment  and  small 
blocks  of  evenly-sampled  flight  data,  indicate  that  unevenly-sampled  flight  data  may  be  analyzed  to  provide 
meaningful  spectral  data. 
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v. 

Figure  12.  (a)  Comparison  between  G' xx(f)  and  Gxx (/)  determined  from  surface  pressure  measurements  recorded 
at  different  times  during  the  ascent,  (b)  Comparison  between  normalized  G'  xx(f)  and  Gxx (/)  versus  normalized 
frequency. 


Figure  13.  (a)  Comparison  between  G' xx(f)  and  Gxx (/)  determined  from  surface  pressure  measurements  recorded 
at  different  times  during  the  ascent,  (b)  Comparison  between  normalized  G' xx(f~)  and  Gxx (/)  versus  normalized 
frequency. 


Table  -  Flight  Conditions  Analyzed 
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Time,  sec 

Freestream 

Mach,  M 

Freestream  Unit  Reynolds 
Number  per  meter 

3.870 

2.14 

4.26E+07 

4.550 

2.61 

4.95E+07 

6.200 

3.43 

5.57E+07 

7.200 

3.25 

4.71E+07 

5.791 

3.41 

5.80E+07 

7.624 

3.19 

4.44E+07 

10.000 

2.90 

3.20E+07 

11.000 

2.83 

2.88E+07 

13.500 

2.67 

2.21E+07 

15.000 

2.58 

1.89E+07 

VI.  Conclusion 

Computing  the  autospectral  density  function  directly  from  an  unevenly  sampled  time  series  results  in  a  distorted 
spectrum.  The  distortion  is  caused  by  nonlinear  sidelobes  that  are  nonzero  only  for  unevenly  sampled  data.  This 
nonlinearity  can  be  removed  by  resampling  the  unevenly  sampled  time  series  onto  an  evenly  spaced  grid  at  the  cost 
of  under  predicting  the  power  of  the  higher  frequency  components  of  the  time  series  data.  A  compensation  method 
was  implemented  to  recover  the  power  lost  by  resampling.  Specifically,  the  compensation  method  developed  in  this 
paper  estimates  the  autospectral  density  function  of  stochastic,  surface  pressure  data  having  red-noise  “like”  spectra 
sampled  unevenly  during  the  HIFiRE-1  flight  experiment.  This  compensation  method  estimates  the  autospectral 
density  function  by  determining  a  correction  factor.  The  coefficients  of  a  first  order  autoregressive  model,  that  holds 
for  unevenly  sampled  data,  are  determined  from  the  measured  time  series  data  using  a  nonlinear  least  squares 
procedure.  The  autoregressive  coefficients  are  used  to  generate  a  time  series  with  the  same  sampling  scheme  as  the 
surface  pressure  data.  This  autoregressive  time  series  and  measured  time  series  are  resampled  onto  an  evenly  spaced 
grid.  The  autospectral  density  function  is  computed  from  the  resampled  autoregressive  time  series,  the  resampled 
measured  data,  and  directly  from  the  autoregressive  coefficients.  A  correction  factor  is  determined  by  taking  the 
ratio  of  the  autospectral  density  function  computed  from  the  autoregressive  coefficients  with  the  autospectral  density 
function  of  the  resampled  autoregressive  time  series.  The  autospectral  density  function  of  the  measured  data  is 
estimated  by  multiplying  the  correction  factor  with  the  autospectral  density  function  of  the  resampled  measured 
pressure  data. 

The  time  series  data  measured  during  the  flight  experiment  has  some  aerothermal  data  that  would  benefit  from 
correlation-based  analysis  such  as  crossflow  instabilities,  rolling  and  pitching  motions  of  the  payload,  and  shock- 
foot  motion  overtop  surface  mounted  transducers.  Future  work  will  investigate  methods  for  preserving  the  phase 
information  of  the  time  series  data  recorded  during  the  HIFiRE-1  flight. 
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Appendix  A:  Least  Squares  Fourier  Transform 

In  this  section,  the  least  squares  periodogram  and  the  least  squares  Fourier  transform  are  developed 
quantitatively  by  summarizing  the  results  given  in  References  17-21.  The  discrete  Fourier  transform  is  defined  as 

FT  {(if)  =  £n=i  yn  exp  (Al) 

where  co  is  the  angular  frequency  given  as  a>=  2nf.  The  classical  periodogram  is  defined  as 

Px(a>)=^\FT(a))\2.  (A2) 

Substituting  Eq.  Al  into  Eq.  A2  gives 

AO)  =  ^l£"=iy„exp(-twtn)|2  (A3) 


or 


AO)  =  ^(En=iynCOS(mtn)]2  +  En=iynSin(wtn)]2}. 


(A4) 
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In  the  derivation  that  follows,  the  time  series  data  y(t )  is  fit  with  the  model  a;-  cos(u>;tn)  +  bj  sin(<u;tn)  using  the 
method  of  least  squares.  The  derived  least  squares  fit  is  then  shown  to  be  equivalent  to  the  classical  periodogram 
defined  in  Eq.  A4.  The  equivalence  includes  the  preservation  of  the  statistical  properties  of  the  classical 
periodogram.  The  least  squares  Fourier  transform  is  then  defined  using  the  least  squares  periodogram. 

The  sum  of  squares  is  given  as 

R  =  In=i{y„  -  [a-j  cos (a>}tn)  +  bj  sin(w;  tn)]}2  (A5) 

For  a  simultaneous  reduction  of  co  =  <o \  oj2,  ...  <op  the  corresponding  parameters  a,  and  bj  are  chosen  such  that  the 
stationary  conditions  =  0.  The  stationary  conditions  are  given  as 

=  EiLibn  COS (cOjtn)  -  dj  cos 2(a)jtn)  -  bj  sin (tuytn)  cos(w;tn)}  =  0  (A6) 

^  =  In=i{yn  Sin(tuytn)  -  Oj  cos ((Ojtn)  sin(o)jtn)  -  bj  sin2 (tuytn)}  =  0  (A7) 


Adopting  the  notation 

CC  =  Xn=iCOS2(a)jtn), 

SS  =  Zn=i  sin2(oJytn), 

cs  =  En=iCos(oj;tn)  sin(ft);-tn), 
YC  =  2n=iynCOs(c0jtn), 
and 

YS  =  Y,n=iyn  sin(ojytn), 

the  stationary  conditions  reduce  to 

r)  R 

—  =  YC  -  d:CC  -  b;CS  =  0 

daj  1  1 

and 

F)  R 

-=YS-  djCS  -  bjSS  =  0. 

dbj  >  > 

Rearranging  Eq.  A13  and  A14  gives 

YC  =  djCC  +  bjCS  =  0, 

YS  =  djCS  +  bjSS  =  0, 


(A8) 

(A9) 

(A10) 

(All) 

(A12) 

(A13) 

(A14) 

(A15) 

(A16) 


YC ]  _ 

CC 

CSl 

'dj' 

-YSi 

-CS 

ss\ 

bJ 

(A17) 


The  two  equations  represented  in  Eq.  A17  are  the  so-called  normal  equations  for  least  squares.  Solving  Eq.  A17  for 
the  dj  and  /?,  coefficients  gives 
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(A18) 


CLj' 

_  rcc 

csr1 

[YC j  _  i 

■  SS 

— CSl 

rci 

Pc 

Lcs 

ss\ 

.ysl  cc-ss-cs2 

.-CS 

CC  i 

YSi 

From  Eq.  A5  and  A18,  the  sum  of  squares  becomes 


A  R  =  R-Z^iyZ 


[YC 

KS] 

■  SS 

— CSl 

rci 

CC-SS- 

-CS2 

.-CS 

CC  i 

-YSi 

(A19) 


At  this  point  AR  is  not  yet  in  the  form  of  the  classical  periodogram.  In  order  to  show  equivalence  between  the  least 
squares  periodogram  and  the  classical  periodogram,  the  model  a;-  cos(u>ytn)  +  bj  sin(u>;tn)  is  replaced  by 
a.j  cos[<Dy(tn  —  t)]  +  bj  sin[<Dy(tn  —  r)]  where  r  is  defined  so  that  CS  =  0  for  each  t„.  Solving 
CS  =  £n=i  cos(u>;tn)  sin(u>;tn)  =  0  for  r  gives 


1^-1  £n=i  sin(2o),-t„) 

—  tan  — s - 7 — - — r 

2  0>j  Xn=iCOs(2a>jtn) 


(A20) 


Equations  A18  and  A19  become 


and 


or 


CLj' 

|  _  r cc  o  r1  \yc-\  _ 

1/cc 

0 

FI 

Pi. 

1 

O 

Co 

Co 

Co 

1 

0 

1/ss. 

AR  =  [YC 


YC z 

AR  =  —  +  —  = 

cc  ss 


By  assuming  that  CC  and  SS  =  N/2,  Eq.  A23  becomes 


AR  * 


(A21) 


(A22) 


(A23) 


(A24) 


AR  «  -  (YC)2  +  —  (YS)2, 

N  V  y  W  V  y 


(A25) 


and 


AR  [£n=iyn  C0s(0Jjtn)]2  +  ^  [Sn=iyn  Sin(oJjtn)]2.  (A26) 

By  comparing  Eqs.  A26  and  A4,  it  is  concluded  that  the  least  squares  periodogram  is  equivalent  to  the  classical 
periodogram. 

While  the  definition  of  r  in  Eq.  A20  conveniently  reduces  the  normal  equations,  establishing  equivalence 
between  the  least  squares  and  classical  periodogram,  it  also  preserves  the  time  invariance  property.  This  can  easily 
be  demonstrated  by  translating  tn  by  T0,  given  as  t„  +  T0.  Under  this  translation,  r  becomes  r  +  T0  and  tu;  (tn  —  r) 
becomes  oij  (tn  +  T0  —  t  —  T0),  which  is  o)j(tn  —  t).  By  maintaining  time  invariance,  the  statistical  properties  of 
the  classical  periodogram  are  identical  to  the  statistical  properties  of  the  least  squares  periodogram. 

With  the  least  squares  periodogram  defined  as 
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(A27) 


^w=®)2+©2. 

it  is  now  advantageous  to  redefine  the  Fourier  transform.  By  comparing  Eqs.  A2  and  A27,  the  magnitudes  of  the 

YC  YS 

real  and  imaginary  components  of  the  Fourier  transform  are  given  as  -=  and  -==,  respectively.  The  sign  of  these 

yjCC  s/SS 

components  is  not  yet  known  since  shifting  each  tn  by  r  introduces  a  phase  shift.  The  phase  shift  can  be  accounted 
for  by  scaling  the  real  and  imaginary  components  of  the  Fourier  transform  by  exp{— ia>j  \tj-  —  r (<"/), )]}  where  tf  =  0 
for  univariate  analysis  and 


.  _  j_y nx  _ Ly  y  t 

lf  -  N„  Ln= 1  lx,n  „  Ln=1  ly,n 


for  bivariate  analysis.  In  practice,  accounts  for  time  series  data  with  different  time  origins. 
To  summarize,  the  least  squares  Fourier  transform  is  given  as 


Y  =  —  ex 
>  V2 


p{-ta)y(t/-T)}(^  +  t|=) 


where  r  is  given  as 


i  -i 

r  =  — tan 

2  (Dj 


a(2a>jtn) 


X^=1cos(2«7tn) 


Substituting  YC,  CC,  YS,  and  SS  into  Eq.  A29  gives 


=  ^expf-toj^^  -  r)}2n 


N  )  ynCQs[ft>j(tK-T)]  ^  yKsin[&)7(tn-T)] 


h;~=1sin2[^(tn-T)] 


(A28) 


(A29) 


(A30) 


(A31) 
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